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Further Refinement of the LEWICE SLD Model 


William B . Wright 
QSS Group, Inc. 

Cleveland, Ohio 44135 

Abstract 

A research project is underway at NASA Glenn Research Center to produce computer software that 
can accurately predict ice growth for any meteorological conditions for any aircraft surface. This report 
will present results from version 3.2 of this software, which is called LEWICE. This version differs from 
previous releases in that it incorporates additional thermal analysis capabilities, a pneumatic boot model, 
interfaces to external computational fluid dynamics (CFD) flow solvers and has an empirical model for 
the supercooled large droplet (SLD) regime. An extensive comparison against the database of ice shapes 
and collection efficiencies that have been generated in the NASA Glenn Icing Research Tunnel (IRT) has 
also been performed. The complete set of data used for this comparison will eventually be available in a 
contractor report. This paper will show the differences in collection efficiency and ice shape between 
LEWICE 3.2 and experimental data. This report will first describe the LEWICE 3.2 SLD model. A semi- 
empirical approach was used to incorporate first order physical effects of large droplet phenomena into 
icing software. Comparisons are then made to every two-dimensional case in the water collection 
database and the ice shape database. Each collection efficiency condition was run using the following four 
assumptions: 1) potential flow, no splashing; 2) potential flow, with splashing; 3) Navier-Stokes, no 
splashing; 4) Navier-Stokes, with splashing. All cases were run with 21 bin drop size distributions and a 
lift correction (angle of attack adjustment). Quantitative comparisons are shown for impingement limit, 
maximum water catch, and total collection efficiency. Due to the large number of ice shape cases, 
comprehensive comparisons were limited to potential flow cases with and without splashing. Quantitative 
comparisons are shown for horn height, horn angle, icing limit, area, and leading edge thickness. The 
results show that the predicted results for both ice shape and water collection are within the accuracy 
limits of the experimental data for the majority of cases. 

Nomenclature 

A p cross-sectional area of particle (nr) 

c chord (m) 

D drag force (kg*m/s 2 ) 

D, e Airfoil leading edge diameter (m) 

d particle diameter (m) 

E total water collection efficiency (dimensionless) 

/ drop frequency (Hz) 

g gravitational constant = 9.8 m/s 2 

h film thickness (m) 

L lifting force (kg*m/s 2 ) 

LWC liquid water content (g/m 2 ) 

m particle mass (kg) 

m water mass flux (kg/m 2 s) 

s surface wrap distance (m) 

s, lower impingement limit by wrap distance (m) 

s ll0 lower limit by wrap distance where (3 = 10 percent (m) 

s u upper impingement limit by wrap distance (m) 
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s uI o upper limit by wrap distance where (3 = 10 percent (m) 

t Time (s) 

V relative droplet velocity (m/s) 

V x x-component of velocity (m/s) 

V y y-component of velocity (m/s) 

x horizontal direction (m) 

x i lower impingement limit by x-location (m) 

x uo lower limit by x-location where (3 = 10 percent (m) 

x„ upper impingement limit by x-location (m) 

Xuio upper limit by x-location where (3 = 10 percent (m) 

x p x-location of water particle (m) 

x first derivative of x-location of water particle with respect to time (x-component of particle 

velocity) (m/s) 

x second derivative of x-location of water particle with respect to time (x-component of particle 

acceleration) (m/s 2 ) 
y vertical direction (m) 

y 0 y-value of the starting locations of collection efficiency trajectories (m) 

y p y-location of water particle (m) 

y first derivative of y-location of water particle with respect to time (y-component of particle 

velocity) (m/s) 

y second derivative of y-location of water particle with respect to time (y-component of particle 

acceleration) (m/s 2 ) 


A. Dimensionless Numbers 


Cl 

c d 

J* 

K 

K l 

k Lii 

K 0 

La 


N, 


f 


Oh 


coefficient of lift = 


A pPa V-/2 
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coefficient of drag = 

v / /2 

dimensionless drop frequency = / ^ 
Mundo splashing parameter = Oh Re"" 


LEWICE splashing parameter = K 
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Re 


Reynolds number = 


We 


pVd 




Weber number = 


p V 2 d 

a 


B. Greek Letters 

a angle of attack (degrees) 

(3 collection efficiency (dimensionless) 

y angle difference between particle velocity vector and airflow velocity vector (radians) 

8 film thickness (dimensionless = hid) 

|1 viscosity (kg/m*s) 

V kinematic viscosity of air (nr/s) 

p density (kg/nk) 

a surface tension (kg/s 2 ) 

9 impact angle (degrees) 

C. Subscripts 

a air 

b bouncing 

/ freezing 

le leading edge 

max maximum value 

n normal direction 

o initial value 

p particle 

5 splashed (or secondary) drop value 

t tangential direction 

term terminal 

tot total 

w water 

x x-dependent 

y y-dependent 

co free-stream property 


I. Introduction 

In 1994, an ATR-72 crashed in Roselawn, Indiana (ref. 1). It has been speculated that accident 
occurred due to the accumulation of ice aft of the deicing boots. Ice formed aft of the boots due to 
impingement of drops greater than 40 pm. Since then, several experimental efforts have been made to 
document supercooled large droplet (SLD) ice shapes and to investigate the underlying physics (refs. 2 
to 7). Based on this experimental work, an empirical model was developed to account for large droplet 
effects in LEWICE (ref. 8). 

This report will address the validation of that model against a database of water collection efficiency 
and ice shape data (refs. 9 to 11). This validation effort mirrors a similar effort undertaken previously for 
the validation of LEWICE for ice shapes (ref. 12). That report quantified the ice accretion prediction 
capabilities of the LEWICE 2.0 software. Several ice geometry features were proposed for comparing ice 
shapes in a quantitative manner. The resulting analysis showed that LEWICE 2.0 compared well to the 
available experimental data. The purpose of this report is to repeat the validation process using the current 
model with the larger databases now available. 
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The report is divided into four sections. The first section will provide a description of the LEWICE 
software and the Navier-Stokes flow solver WIND (ref. 13) that was used for this effort. The second 
section will provide a description of the LEWICE collection efficiency model with emphasis on breakup 
and impact physics. It will also describe the modified equations including analysis and observations from 
tests performed in the NASA icing research tunnel (IRT). The third section will describe the experimental 
data and the parameters used for quantifying the comparisons. The last section will provide validation 
results along with a statistical comparison of those parameters with the available experimental data. 

II. Computational Tools 

A. LEWICE 

The computer program, LEWICE, embodies an analytical ice accretion model that evaluates the 
thermodynamics of the freezing process that occur when supercooled droplets impinge on a body. The 
atmospheric parameters of temperature, pressure, and velocity, and the meteorological parameters of 
liquid water content (LWC), droplet diameter, and relative humidity are specified and used to determine 
the shape of the ice accretion. The surface of the clean (un-iced) geometry is defined by segments joining 
a set of discrete body coordinates. The software consists of four major modules. They are 1) the flow field 
calculation, 2) the particle trajectory and impingement calculation, 3) the thermodynamic and ice growth 
calculation, and 4) the modification of the current geometry by addition of the ice growth. 

LEWICE applies a time-stepping procedure to “grow” the ice accretion. Initially, the flow field and 
droplet impingement characteristics are determined for the clean geometry. The ice growth rate on each 
segment defining the surface is then determined by applying the thermodynamic model. When a time 
increment is specified, this growth rate can be transformed into an ice thickness and the body coordinates 
are adjusted to account for the accreted ice. This procedure is repeated, beginning with the calculation of 
the flow field about the iced geometry, then continued until the desired icing time has been reached. The 
results shown in this report are from version 3.2 of LEWICE (ref. 14). 

B. WIND 

WIND is a structured, multi-zone, compressible flow solver with flexible chemistry and turbulence 
models. Zonal interfaces may be abutting or overlapped, allowing the flexibility to treat complex systems 
moving relative to one another. WIND is a computational platform that may be used to numerically solve 
various sets of equations governing physical phenomena. Currently, the software supports the solution of 
the three-dimensional Euler and Navier-Stokes equations of fluid mechanics, along with supporting 
equation sets governing turbulent and chemically-reacting flows. 

WIND is a product of the NPARC Alliance, a partnership between the NASA Glenn Research Center 
(GRC) and the Arnold Engineering Development Center (AEDC) dedicated to the establishment of a 
national, applications-oriented flow simulation capability. The Boeing Company has also been closely 
associated with the Alliance since its inception, and represents the interests of the NPARC User’s 
Association. 

C. ICEG2D 

ICEG2D (ref. 15) is an automated grid generation program and scripting interface. The grid 
generation capability contains algorithms for producing single block “C” grids. The scripting portion of 
the software provides an interface between LEWICE and WIND. Flow field information from WIND is 
processed through the WIND utility GMAN and formatted for use in LEWICE. LEWICE then sends the 
completed ice shape back to ICEG2D for the next icing time step. 
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III. Collection Efficiency Physics 

Collection efficiencies are calculated in LEWICE through the use of a particle trajectory analysis. 
Droplets are released from a point in the freestream flow and tracked through the flow field using the 
following equations: 

mx p = —D cos y - L sin y my p = -D sin y + L cos y - mg 

y„-V v of 2 o V 2 i 7 7 

where Y = a tan . , D = c d ^—A p , L = c l ^ L —A p , V = ^(x p -V x ) +(y p ~V y ) 2 

X p V x L 2. ' 


The initial release point in the x-direction is determined by finding an x-location where all of the 
velocities in a vertical sweep are within 0.1 percent of the freestream value. The initial release point in the 
y-direction is determined from the angle of attack. The initial velocity is assumed to be the terminal 
droplet velocity, which is given by 


c d Re 


2 

term 


4 gd (Pw~P q) 


, _ Vfgfm d 

where R e term = 

^ a 

Each droplet is then tracked until it either hits the airfoil or reaches the trailing edge. After the first 
trajectory ends, the next particle is released from a higher or lower starting point in an attempt to hit the 
surface. This process is continued until there exists at least one drop that passes above the airfoil and one 
that passes below. 

Impingement limits are then found using a standard bisection search algorithm. Using the coarse 
limits found in the prior step, LEWICE starts a drop halfway between these limits. Based upon the end 
result of that trajectory, the next drop is released halfway between the current starting point and the 
starting point of either the upper or lower coarse limit. The coarse limit is then refined based upon the 
trajectory results. This process is repeated until the starting point of a drop that hits the airfoil and the 
starting location of a drop that misses the airfoil is within 10“ 5 (dimensionless) The bisection is then 
repeated to find the second impingement limit. 

Collection efficiency is then determined by sending a user-determined number of trajectories 
uniformly spaced between the impingement limit starting locations. The starting and ending locations of 
these trajectories are stored. Collection efficiency is then calculated by the following definition: 

n dy 0 
ds 


A. Literature review 

A previous report published a literature review of large droplet physics (ref. 8). In that report, the 
additional physics that were found to be of most importance were drop breakup and drop splashing. While 
many of the other physical effects can be studied, this report will focus on those that have the most impact 
on collection efficiency. 
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1. Drop Breakup 

If a large drop moves at a high enough velocity, it can break up due to shear. Breakup occurs when 
the drop passes a critical Weber number. Values for this critical Weber number vary widely in the 
literature. 

The Weber number is given by: 


We p 


P gV 2 d 

CT 


For water drops falling at their terminal velocity, the critical Weber number (based on air density) is 
approximately 10. For water drops accelerated by a shock wave, a value of 6.5 is given. Krzeczkowski 
(ref. 19) and Hsiang (ref. 20) each measured droplet breakup for shear-induced flows and reported values 
ranging from 10 to 20. Ibrahim et al. (ref. 21) provided a more detailed analysis of the droplet 
deformation and breakup using a Taylor analogy model. The Weber number of each trajectory was output 
from LEWICE for the case described above to investigate this effect. 

An empirical relationship found in Hsiang and Faeth (ref. 20) was added to LEWICE in order to 
estimate the reduction of collection efficiency due to breakup. In their model, droplets will start to break 
up when the critical Weber number is greater than 13. This Weber number is based on the air density as 
defined earlier in this report. Since droplet breakup occurs rapidly compared to the trajectory time step, 
breakup is considered to be instantaneous. Dai and Faeth (ref. 22) produced some excellent photographs 
of the breakup process using pulsed shadowgraphy and holography. 

Secondary particle size is given by the following equation: 


d s = 6.2| 


( p V4 

r iv 

l.p« . 


R e>v 



where Re M , = 

Fw 

This correlation can be applied to an Eulerian system as well as the Lagrangian tracking system used 
by LEWICE. However, it would be necessary in an Eulerian system to solve coupled sets of equations for 
each drop size generated. In the Lagrangian system, the smaller drop size is simply tracked from the 
breakup location. An empirical relationship was chosen to assess the importance of breakup to the 
collection efficiency. If breakup is not important then there is no need to implement the more complicated 
droplet deformation and breakup (DDB) model described by Ibrahim (ref. 21). 

2. Drop Splashing 

The literature search performed for this report confirms that the primary assumption used by 
LEWICE 2.0 that was invalidated for SLD was the assumption that all drops that strike the surface 
impinge, thus neglecting splashing and/or bouncing of drops. One of the earliest detailed experimental 
studies was performed by Stow and Hadfield (ref. 23) who reported on the impact of water drops on a dry 
surface. Macklin and Metaxas (ref. 24) reported a similar study that also used ethanol and glycerol to 
study the effect of different fluid properties. Jayarante and Mason (ref. 25) looked at bouncing and 
splashing of raindrops impinging at various angles on dry surfaces and films. Wright (ref. 26) developed a 
theoretical splash model for raindrops for the purpose of modeling soil erosion. 

Harlow and Shannon (ref. 27) solved the Navier-Stokes equations for the impact of a single drop on a 
dry surface or film. A more recent work was performed by Yarin and Weiss (refs. 28 and 29) who 
proposed a splashing model as a type of kinematic discontinuity. Other works include Rein (ref. 30) who 
provides a review of several papers, including phenomena such as bouncing along with splashing and 
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coalescence, and Chandra and Avedisian (ref. 31) who documented photographically the droplet structure 
during the deformation process. 

Computationally, a detailed physical model of droplet splashing would require solving the Navier- 
Stokes equations for each droplet impact using a Volume of Fluid (VOF) model such as that described by 
Hirt (ref. 32). Examples of this type of calculation were reported by Trapeaga and Szekely (ref. 33) as 
well as Tan and Papadakis (ref. 2). Current computational capabilities usually limit this approach to single 
drop calculations. In a typical icing encounter, thousands of droplet impacts are recorded per second, 
making this type of analysis prohibitively expensive. An empirical or semi-empirical approach is 
therefore necessary. 

Experimental studies by Mundo, Sommerfeld, and Tropea (refs. 34 to 36) examined droplet-wall 
collisions and correlate splashing in terms of Reynolds number and Ohnesorge number: 


Oh = 


VWe 

Re 





The Reynolds and Ohnesorge numbers are based on the liquid (water) properties and the component 
of the impact velocity normal to the surface. Based on the results of their experiment, splashing occurs if 
the factor K = Oh*Re' 25 is greater than 57.7. 

The Mundo papers also provide a characterization of the size, velocity, and direction of the splashed 
particles. Their later references provide an empirical splashing model that can be used in Lagrangian 
tracking schemes. The empirical formulas calculate splashed drop size, splash velocity, splash angle, and 
deposited mass as functions of the incoming parameters. In last year's report, several problems were noted 
with the empirical relationships developed at low speeds (<30 m/s). The following section provides a 
summary of the current splashing model. 

3. Summary of LEWICE Splashing Model 

The parameters needed for an empirical splashing model are the splashing threshold, the splashed 
drop size (or drop size distribution), the splashed velocity (or a distribution of splash velocities), the 
splash angle (or a distribution of splash angles), and the amount of splashed mass. The number of 
splashed particles is needed only if all splashed particles are to be tracked. The models presented in the 
literature typically provide these variables as a ratio to the incoming parameters. 

The droplet experiments mentioned previously were developed using similar experimental techniques 
and therefore had similar ranges of applicability. The upper limit on drop size and velocity were 340 pm 
and 30 m/s respectively. Droplet frequency and film thicknesses were in the range expected for icing 
encounters except at the lower range. Splash data exists for dry surfaces (5 = 0) and for film thicknesses 
of 0.3 < 8 < 3. The applicability of these models to very thin films is unknown. Similarly, droplet 
frequencies are lower in the SLD range due to the lower water contents as well as the higher drop sizes. 
However, the correlations are well behaved in these limits and tend toward limiting values. 

Previous reports used a modified version of empirical relationships reported by Trujillo et al. (ref. 

38). In addition, a bouncing model was incorporated into LEWICE based on physical arguments that 
since droplet impact phenomena at the high velocities encountered in aircraft icing have not previously 
been studied in the literature, it is reasonable to assume that additional phenomena could be present at 
these velocities. In the current model, both effects were combined into a single correlation. The 
correlation was determined by matching the collection efficiency results to the experimental collection 
efficiency database. All 114 cases in this database were used to produce this calibration. This process was 
made possible through the use of a Visual Basic script that ran the cases and made collection efficiency 
charts in Excel for all 1 14 runs. The procedure used was to run all the cases in batch mode and then 
compare the collection efficiencies for all of the runs. The splashing correlation was then adjusted to 
match the largest number of cases possible. Grid-independent solutions of the Navier-Stokes equation 
were generated using WIND and then used for the flow solution inputs. Because the airfoil shape does not 
change during the impingement tests, the flow solutions need only be generated once. Since the collection 
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efficiency database was run at one airspeed (175 mph), several runs were also performed at higher speeds 
and at larger drop sizes in order to assess trends. The correlations that best matched the collection 
efficiency database can be expressed by: 

Splashing threshold: K Ln > 200 


where 


K 


L,n 



^ s _ g rj 2 g— 0 . 028 \K 


0.05 < — < 1; ^ = O.7[l-sin0 o l [l-exp(-0.0092*(tf in -200)) 

d n m n v ’■ 


V, V 

-^- = 1.075-0.00250 ' — = 0.3 - 0.0029„ 

V O’ y O 

’ tfi r n 


It should be noted that if droplets impinge perpendicular to the surface (9o = 90°) then no splashing can 
occur. While this result is the opposite trend of the low velocity experiments, a large number of SLD 
cases show very little mass loss at the leading edge. Mass loss in the collection efficiency data occurs 
primarily at the low impact angles near the impingement limits except for the MVD = 92 pm and to a 
lesser extent the MVD = 79 pm data. A feature was added to LEWICE to track the trajectories of the 
splashed particles and the trajectories of particles after breakup. This process was described in a report by 
Rutkowski et al. (ref. 7). Mass losses due to splashing and subsequent re-impingement are also calculated. 
Since it is necessary to extrapolate from the experimental data for icing encounters, care must be taken, 
especially in high velocity impacts. The equations provided above indicate that droplet velocity (which is 
related to airspeed) is at least as important to the splashing process as drop size. The next section will 
describe the experimental data in more detail. 


IV. Experimental Data 


A. Collection Efficiency Data 

The experimental data used for this comparison was taken from several tests performed in the IRT 
(refs. 9 to 11). The tests were performed by personnel from Wichita State University, NASA, and Boeing 
using a spray system designed for short duration sprays. The water spray contains a known concentration 
of blue dye and the models are covered with a heavy weight blotter paper. The amount of dye was then 
measured via reflectance spectroscopy using a CCD camera. Collection efficiencies were determined on 
1 1 different clean airfoils and 7 different iced airfoils. Drop size distributions from MVDs of 1 1 to 
236 pm were used and angle of attack was varied from 0° to 8°. 

The database consists of a wide variety of airfoils that were selected to provide a broad range of 
airfoil types for validation purposes. A previous report (ref. 8) contained a summary of the conditions and 
airfoils used in the collection efficiency tests. The drop size distributions were measured during the 
experiment and were modeled in LEWICE using 27-bin distributions. The estimated experimental 
repeatability for this data ranged from 10 to 30 percent for the older data sets. An uncertainty analysis 
presented in the reports on the experimental data showed that uncertainty in spray time, dye 
concentration, spray pressures, tunnel velocity, and cloud unsteadiness due to turbulence accounted for a 
14 percent variation in maximum collection efficiency values. In the following comparison with predicted 
results, the 14 percent value was used as representative of the experimental error for all parameters. 
Specific variability of individual parameters could not be defined (ref. 39). 
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B. Ice Shape Data 

The experimental data described in this paper are the result of a wide variety of tests performed in the 
NASA Glenn Icing Research Tunnel (IRT) in recent years (refs. 17 to 24). Eight airfoils and a cylinder 
were selected for this comparison. These airfoils and the accompanying ice shapes represent the complete 
set of publicly available data that has been generated in the IRT and digitized for single element airfoils. 
There are a total of 2146 IRT runs analyzed for this validation report, of which 460 are repeats of 
previous runs in the IRT. There are 865 digitized tracings at off-centerline locations for a total of 301 1 
experimental ice shapes. This value is nearly triple the number of ice shapes from the previous validation. 
There are 741 cases with an MVD greater than 40 pm. Of these cases, 111 are rime or roughness (short 
duration) conditions and 59 cases use SLD conditions for only part of the run. The eight airfoils are as 
follows: a modified NACA23014, a Large Transport Horizontal Stabilizer (LTHS), a GLC305 business 
jet airfoil, a NACA0012, a modified NACA4415, a NACA 0015, a NACA23012, and a NLF00414 
laminar flow general aviation wing. These models were described in an earlier report on the LEWICE 2.0 
validation effort (ref. 12). 

The data is taken in the IRT by cutting out a small section of the ice growth and tracing the contour of 
the ice shape onto a cardboard template with a pencil. The pencil tracing is then transformed into digital 
coordinates using a stylus on a digitizing table. A flat-bed scanner with digitizing software has been 
available to accelerate the data acquisition process. For any given IRT test run, up to five span wise 
sections of the ice shape are traced and digitized in this manner. There are several steps within this 
process that can potentially cause experimental error. Those that can be quantified by the current 
technique are the span wise variability, the repeatability error, and errors involved in the tracing technique. 


V. Results and Comparison Methodology 

A. Collection Efficiency Results 

The comparison process for collection efficiency mimics the process used in a previous report 
(ref. 12) for comparing ice shapes. The following six parameters were chosen to characterize the 
collection efficiency results: maximum collection efficiency ((3 max ), total collection efficiency (E) which is 
proportional to total water catch rate, upper and lower impingement limits (y, and y) and the location on 
the upper and lower surface where collection efficiency reached 10 percent (y,io and yi 0 ). These 
parameters are shown on a representative collection efficiency curve in figure 1 below. The parameters 
are labeled on the figure for both the experimental and numerical results. The example shown is for the 
GLC305 airfoil at 1° AOA and 168 pm MVD. While this figure depicts the impingement distance based 
on wrap distance, the impingement limit comparisons were done using the x-distance from the leading 
edge. For the ice shape cases, x/c = 0 represents the leading edge of the clean airfoil geometry. 

The six parameters were chosen as being of interest in designing deicing systems. The maximum 
collection efficiency and total water catch are critical parameters for evaporative systems. Impingement 
limits are used extensively for designing the limits of ice protection systems. The 10 percent impingement 
limits were included in this parameter list since some use this location to define the limits of ice 
protection. This assumption may be due to the insensitivity of an airfoil to small amounts of ice or based 
on user assumptions of the software's accuracy. The six parameters were measured for the experimental 
data and for each of the predicted runs. The values were then compared by absolute value, percent 
difference and for impingement limits by percent chord difference. Results will be shown for the average 
comparisons as well as a statistical analysis of the software accuracy. 

Comparisons of the six parameters were made between LEWICE and the experimental data. Since the 
current correlation is based upon the experimental database it is compared against, this report will only 
show a limited number of results in order to concentrate on the ice shape comparison. Since the splashing 
and breakup model is not complete, it can be activated or deactivated by the user. Although all effects are 
modeled, this report will refer to the splashing, breakup, and bouncing model as simply the splashing 
model since droplet splashing had a much larger effect on the results than droplet breakup. LEWICE was 
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ran four times for each condition in the database. In the first option, LEWICE was run without splashing 
effects and using Navier-Stokes flow solutions from WIND. This option was labeled “NS-NS” (Navier- 
Stokes, no splashing). In the second option, splashing was considered along with the Navier-Stokes flow 
solver. This option was labeled “NS-S”. The third option used the embedded potential flow solver and 
ignoring splashing and breakup. This option is labeled “PF-NS” (potential flow, no splashing) in some 
charts. In the fourth option, LEWICE was run using the potential flow solver and considering splashing. 
This option is labeled “PF-S”. In the third option, the Navier-Stokes flow solver WIND replaced the 
potential flow solver. All cases were run with 27 bin drop size distributions as measured in the 
experiment. 

1. Total Mass Flux 

The total mass flux was obtained by integrating the collection efficiency values with respect to the 
LWC *V r 

wrap distance, m = * J fids . Percentage differences were calculated and tabulated with the other 

results. Statistical averages were also measured for this variable. Figure 2 shows the prediction of mass 
flux for the four LEWICE scenarios. The solid bars show the average difference between the prediction 
and experiment while the vertical error bars show the variability (standard deviation) of the results. 

This figure best demonstrates the improvement achieved by using the splashing model. Changing the 
flow solver to Navier-Stokes had only a small effect on water catch while there is a significant over- 
prediction of water catch if splashing is not included. The current splashing model corrects the problem 
shown previously of removing too much mass in the lowest drop size range. The model shows a 
consistent improvement for SLD drop sizes of nearly four-fold improvement in accuracy for water mass 
flux (from a difference of 60 or 80 to a 20 percent difference). The average difference in the SLD range is 
now comparable to the accuracy for the appendix C cases. Since the current correlation is based upon the 
same data set that it is compared against, a good comparison is assured for all drop sizes. During the 
calibration process, correlations were found that had a lower average error for the parameters selected, but 
the amount of under prediction was deemed too great. It was desired to have a correlation that was 
slightly conservative when an exact match to the experimental data could not be found. The next section 
will show how this calibrated model performs against the ice shape database. 

B. Ice Shape Results 

This section describes the methodology used to make the quantitative measurements on experimental 
and predicted ice shapes. This methodology has been incorporated into a software utility called THICK 
that calculates and outputs the parameters described. This software was created in order to process the 
large number of ice shapes presented in this report. This program reads two geometry files: one for the 
clean airfoil and one containing an ice shape. This utility has been documented in the LEWICE 3.0 User 
Manual (ref. 2). This software has recently been revised to more accurately capture ice features of 
interest. In particular, ice shapes are first segregated into the following categories: roughness shapes (low 
exposure times), rime shapes, mixed shapes (0.6 <Nf< 1), and glaze (Nf <0.6). Roughness shapes use 
icing limit for comparison. Rime shapes use area, maximum thickness, and icing limit. The glaze and 
mixed shapes are categorized using icing limit, area, horn height, horn thickness, and leading edge 
thickness. The parameters were defined in the previous validation report, although area is now defined as 
the true area rather than the integrated thickness. The parameters are diagrammed in the figure below. 

A total of 3049 LEWICE cases were performed in this validation effort. There were 1626 unique 
cases with one spray condition and 59 for cases with multiple spray conditions. The 944 appendix C cases 
were all run with a monodispersed drop size. The 682 SLD cases were run three times. The first set used a 
monodispersed drop size equal to the MVD. The second set of SLD runs used a 27 bin drop size 
distribution. Since drop size distributions in the tunnel are only known at calibration points, the drop size 
distributions for these runs were estimated from the known distributions. The normalized distribution 
used for these runs is shown in figure 4. This normalized distribution is then multiplied by the MVD to 
obtain the droplet distribution used for the LEWICE runs. The last set of runs used the 27 bin drop size 
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distributions with the splashing effects turned on. Since splashing effects depend so greatly on the drop 
size, it was not considered appropriate to use the splashing effects with a monodispersed drop size. 

Due to the large number of cases, results were obtained using automated scripts. For the LEWICE 
cases, runs were produced using a utility program called LewCon. This program was created for this 
validation effort. It was written by Aaron Armstrong from Brigham Young University at Idaho and later 
modified by this author. The utility is included with the latest LEWICE release. The program loads 
conditions from an Excel spreadsheet and uses those values to create the LEWICE input files. The 
program also runs LEWICE and places the results, including charts, into individual Excel files for each 
case run. The experimental data was added to these charts using several Excel macros developed by this 
author and by Timothy Mosig from the Worchester Polytechnic Institute. The Excel macros also ran the 
THICK utility on the LEWICE shapes and the experimental ice shapes to generate the ice shape 
measurements presented in this report. The following sections will describe the results of that process by 
showing comparisons of different ice shapes with LEWICE. 

C. Icing Limits 

The icing limits are the chordwise locations on the ice shape on the upper and lower surface where 
the ice shape merges with the airfoil. Both the wrap distance from the leading edge and the .r-di stance are 
recorded for each icing limit. The results presented here are for the wrap distance values. 

Figure 5 shows the results of these measurements for both the experimental ice shapes and for 
LEWICE. In this comparison, each measurement on an experimental shape is compared against the 
average measurement for that condition. It does not matter whether the experimental shape is a repeat run 
or a tracing at a different location since the previous validation report showed that the variation is very 
similar. The LEWICE comparison is the difference between the measurement on the LEWICE shape and 
the experimental average value. If only one experimental ice shape was available for a given condition, 
the measurement on that experimental ice shape was taken as the average value for the condition and 
compared to LEWICE. These results are presented as a percentage of chord in order to normalize the 
results for different cases. This figure shows that the experimental variation in the lower icing limit is 2 
percent of chord while the LEWICE result lies within 4.5 percent of chord from the experimental average 
value. This result uses the absolute error for each case in order to compute the average and represents an 
improvement over the previous validation effort. Contrary to popular belief, in the majority of cases 
LEWICE underpredicts rather than overpredicts the icing limit as compared to the experimental data. This 
result can likely be attributed to the use of a monodispersed drop size when obtaining the predicted result. 
Additionally, there were no discernable trends in the data due to ice type (rime, mixed, glaze). Figure 6 
shows the comparison just for the SLD database. Use of a droplet distribution, while physically more 
accurate, makes the agreement with experiment slightly worse. 

D. Maximum Ice Thickness (Horns) 

The current release of the THICK utility has an improved methodology for finding the glaze ice 
horns. The routine first calculates the ice thickness normal to the surface at each location. Then the 
program finds the “convex hull” of the ice shape and selects the two largest thicknesses on the convex 
hull as the ice horns. The convex hull are those points on a geometry that create a convex surface. In order 
to compare different conditions with different chord lengths and accretion conditions, the individual ice 
thicknesses were non-dimensionalized by the maximum accumulation thickness as given in the equation 
below. 


(. LWC)(v)(Time ) 

P ice 


Figure 7 shows the dimensionless difference in ice thickness for the three ice thickness measurements 
made in this report. The previous validation showed that repeatability differences and spanwise variations 
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were on the same magnitude, so they were not plotted separately in this report. This figure shows that the 
max. thicknesses can be measured to within 18 percent experimentally and that the average difference for 
the LEWICE cases is 

20 percent for max. thickness. Both the experimental valuation and the LEWICE predicted values show a 
higher difference for horn height than those reported in the previous validation effort. This occurred 
solely due to errors in the current automated procedure. The large error bars for the mixed ice shape 
comparison are an indication of this effect. The utility program THICK does not always extract the 
correct location for the ice horns. However, in order to correct this result and present a proper comparison 
would require manual assessment of each horn selection. This process could not be completed in time for 
this report. Figure 8 shows the comparison for just the SLD data. Once again although the use of a droplet 
distribution and a splashing model are considered more physically accurate, that does not necessarily 
translate to actual comparisons. 

E. Ice Area 

The iced area is calculated by finding the area of the iced surface through integration, then subtracting 
the area of the clean airfoil. A unified comparison of ice area for cases with different absolute areas also 
poses a problem. In this report, the area difference has been nondimensionalized by the maximum 
accumulation thickness given earlier and by the airfoil thickness. It should be noted that the absolute 
values for ice area are maintained in the Microsoft® Excel (Redmond, WA) spreadsheet so that the users 
of this data can make their own comparisons. 

Figure 9 shows the results for the ice area comparison for all cases and for SLD conditions. This 
figure shows that the experimental difference in ice area is less than 10 percent while for the LEWICE 
results the variation is approximately 16 percent. Once again there are no discernable trends in the 
comparison based on ice type. 

F. Angle at Max. Thickness (Horn Angle) 

The horn angle was defined as the angle between two line segments. The first line segment is a 
horizontal line which goes through the leading edge of the airfoil (minimum x-location at AOA = 0°). The 
second line segment is defined by the leading edge of the clean airfoil and the ice horn location. This 
angle was measured for all ice shapes whether or not they fit the classical definition of having a glaze ice 
horn. Because of this, the ice shapes were first separated into four categories: roughness, rime, mixed, and 
glaze shapes. The separation of the data into different categoriesallows for a truer comparison of results. 

Figure 10 shows the variation in max. thickness angle for LEWICE and for the experimental 
categories described earlier. Results are presented in degrees. This figure shows that the variation in the 
experimental data is 26° for the upper angle, and 31° for the lower angle. The LEWICE difference from 
the experimental average are essentially the same as the experimental variation. Once again, however, this 
valuation is primarily due to errors in the automated process and not due to actual differences in shapes. 
Qualitative assessment indicates that the relative comparison is accurate. That is, the difference in 
prediction from LEWICE to the experimental average is on the same order as the variability in the data 
itself. Figure 11 shows the comparison with just the SLD data. Similar to the other measurements, the 
LEWICE SLD results are slightly better than the overall comparison and is made slightly worse by 
introducing a distribution with splashing. This could also indicate that the splashing model could be better 
calibrated for ice shapes as well as for water impingement. However due to the large number of ice shape, 
this was not practical. 

G. Overall Assessment 

Once the individual measurements are taken for each ice shape, it becomes useful to create an overall 
assessment of the ice shape prediction. Since each measurement is different, several methods could be 
used to assess the overall difference between two ice shapes. Eight of the 10 measured values presented in 
this report have been nondimensionalized. Angles do not have a characteristic measure to use for 
nondimensionalization, so the three angle criteria are reported in degrees. The overall assessment was 
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determined by an average of the eight individual dimensionless values and the three angle criteria in 
degrees. 

Figure 12 shows the comparison of the overall assessment for the experimental variability and for 
LEWICE. This calculation shows an average overall difference of approximately 12 percent for the 
experimental database and 14 percent for LEWICE. The standard deviation is 16 percent for both the 
experimental data and for the LEWICE results. In order to determine if this simple average is a good 
assessment of the variation, plots were made of the average variation for the experimental shapes and for 
the LEWICE shapes. 

Figure 13 shows an example of two ice shapes which are near the overall experimental average. This 
plot shows the spanwise variability from a data point in the LTHS database. The qualitative comparison 
of these two ice shapes suggests that the overall assessment parameter is a reasonable approximation. The 
comparison also shows the problem with automatically generating ice measurement parameters. Although 
the two shapes are very simliar, THICK chooses locations for both the upper and lower horn for these 
shapes that are very different which affects the quantitative analysis. Similarly, figure 14 shows an 
example which is at the average variation for the LEWICE cases. The qualitative assessment of this 
comparison also agrees with the overall assessment parameter used. 

VI. Conclusions 

A review of the physical effects of large droplet phenomena confirmed that droplet splashing is the 
primary physics that is not modeled in icing software in this drop size regime. Empirical models for these 
phenomena were presented. The resulting effects of the additional physics on water collection in 
LEWICE were presented. 

The results show that the empirical splashing models in the open literature only accounted for 
splashing effects near the leading edge. The models, all of which are based on the velocity component 
normal to the surface, could not explain the measured decrease in water collection near the impingement 
limits. Additionally, it was necessary to extrapolate values from that empirical data for this effort. It is 
hypothesized that droplets that impact at high total velocities (> 30 m/s) but low impact angles (< 20°) 
may rebound in a manner different than the current experimental data at the lower velocities. A new 
model based on physical arguments was proposed to account for the additional physical phenomena that 
occur in the SLD regime. 

The model was calibrated against a database of collection efficiencies. A comprehensive comparison 
was performed that compared maximum collection efficiency, total water mass flux, upper and lower 
impingement limits, and upper and lower 10 percent limits. Comparisons were made with and without the 
new model to demonstrate its effectiveness. Comparisons were also made to assess the effect of using a 
Navier-Stokes flow solver and to assess the effect of other adjustments including the use of 27-bin drop 
size distributions instead of 10-bin distributions. 

The accuracy of the model was demonstrated against a voluminous database of ice shapes measured 
in the IRT. Over 3000 experimental ice shapes were used in this comparison. Ice shape measurements 
were made for upper and lower icing limit, upper and lower horn height, leading edge minimum 
thickness, area and upper and lower horn angle. Due to the large number of cases, ice shape generation 
and comparison were performed using several scripts which limited the amount of time that could be 
spent manually evaluating the comparisons. The measurements showed that LEWICE predictions 
compared to the experimental averages were on the same level of accuracy as the tunnel variability for all 
parameters, even in the SLD regime. The overall comparison showed the LEWICE results at 16 percent 
with a standard deviation of 16 percent while the experimental variation was 16 percent with a standard 
deviation of 16 percent. 
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Figure 1. — Sample collection efficiency curve. 
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Figure 2. — Comparison of total mass flux. 
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Figure 3. — Ice shape measurements. 
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Figure 4. — Drop size distribution for SLD cases. 
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Figure 5. — Icing limit comparison. 
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Figure 6. — SLD icing limit comparison. 
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Figure 7. — Horn height comparison. 
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Figure 8. — SLD horn height comparison. 
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Figure 10. — Horn angle comparisons. 
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Figure 11. — SLD hom angle comparison. 
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Figure 12. — Overall assessment. 
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Figure 14. — Average LEWICE variation. 
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